BCR Work-Precision Diagrams
Samuel Isaacson and Chris Rackauckas
The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff chemical reaction network modeling the BCR signaling network from Barua et al.. We use ReactionNetworkImporters to load the BioNetGen model files as a Catalyst model, and then use ModelingToolkit to convert the Catalyst network model to ODEs.
using DiffEqBase, OrdinaryDiffEq, Catalyst, ReactionNetworkImporters, Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq, LSODA, TimerOutputs, LinearAlgebra, ModelingToolkit, BenchmarkTools gr() datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") const to = TimerOutput() tf = 100000.0 # generate ModelingToolkit ODEs @timeit to "Parse Network" prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) rn = prnbng.rn @timeit to "Create ODESys" osys = convert(ODESystem, rn) tspan = (0.,tf) @timeit to "ODEProb No Jac" oprob = ODEProblem(osys, Float64[], tspan, Float64[]) @timeit to "ODEProb DenseJac" densejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true)
Parsing parameters...done
Creating parameters...done
Parsing species...done
Creating species...done
Creating species and parameters for evaluating expressions...done
Parsing and adding reactions...done
Parsing groups...done
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100000.0)
u0: 1122-element Vector{Float64}:
299717.8348854
47149.15480798
46979.01102231
290771.2428252
299980.7396749
300000.0
141.3151575495
0.1256496403614
0.4048783555301
140.8052338618
⋮
1.005585387399e-24
6.724953378237e-17
3.395560698281e-16
1.787990228838e-5
8.761844379939e-13
0.0002517949074779
0.0005539124513976
2.281251822741e-14
1.78232055967e-8
@timeit to "ODEProb SparseJac" sparsejacprob = ODEProblem(osys, Float64[], tspan, Float64[], jac=true, sparse=true) show(to)
──────────────────────────────────────────────────────────────────────────
────
Time Allocations
─────────────────────── ────────────────────
────
Tot / % measured: 1.36h / 99.9% 1.81TiB / 100.0%
Section ncalls time %tot avg alloc %tot
avg
──────────────────────────────────────────────────────────────────────────
────
ODEProb DenseJac 1 1.20h 88.1% 1.20h 1.60TiB 88.4% 1.6
0TiB
ODEProb SparseJac 1 505s 10.3% 505s 200GiB 10.8% 20
0GiB
ODEProb No Jac 1 42.0s 0.9% 42.0s 7.25GiB 0.4% 7.2
5GiB
Create ODESys 1 19.1s 0.4% 19.1s 4.05GiB 0.2% 4.0
5GiB
Parse Network 1 15.6s 0.3% 15.6s 2.52GiB 0.1% 2.5
2GiB
──────────────────────────────────────────────────────────────────────────
────
@show numspecies(rn) # Number of ODEs @show numreactions(rn) # Apprx. number of terms in the ODE @show length(parameters(rn)) # Number of Parameters
numspecies(rn) = 1122 numreactions(rn) = 24388 length(parameters(rn)) = 128 128
Time ODE derivative function compilation
As compiling the ODE derivative functions has in the past taken longer than running a simulation, we first force compilation by evaluating these functions one time.
u = ModelingToolkit.varmap_to_vars(nothing, species(rn); defaults=ModelingToolkit.defaults(rn)) du = copy(u) p = ModelingToolkit.varmap_to_vars(nothing, parameters(rn); defaults=ModelingToolkit.defaults(rn)) @timeit to "ODE rhs Eval1" oprob.f(du,u,p,0.) @timeit to "ODE rhs Eval2" oprob.f(du,u,p,0.) densejacprob.f(du,u,p,0.) sparsejacprob.f(du,u,p,0.)
We also time the ODE rhs function with BenchmarkTools as it is more accurate given how fast evaluating f is:
@btime oprob.f($du,$u,$p,0.)
37.069 μs (3 allocations: 576 bytes)
Now we time the Jacobian functions, including compilation time in the first evaluations
J = zeros(length(u),length(u)) @timeit to "DenseJac Eval1" densejacprob.f.jac(J,u,p,0.) @timeit to "DenseJac Eval2" densejacprob.f.jac(J,u,p,0.)
ERROR: syntax: expression too large
Js = similar(sparsejacprob.f.jac_prototype) @timeit to "SparseJac Eval1" sparsejacprob.f.jac(Js,u,p,0.) @timeit to "SparseJac Eval2" sparsejacprob.f.jac(Js,u,p,0.) show(to)
──────────────────────────────────────────────────────────────────────────
────
Time Allocations
─────────────────────── ────────────────────
────
Tot / % measured: 1.57h / 99.6% 1.82TiB / 99.9%
Section ncalls time %tot avg alloc %tot
avg
──────────────────────────────────────────────────────────────────────────
────
ODEProb DenseJac 1 1.20h 76.5% 1.20h 1.60TiB 87.8% 1.6
0TiB
SparseJac Eval1 1 700s 12.4% 700s 6.21GiB 0.3% 6.2
1GiB
ODEProb SparseJac 1 505s 9.0% 505s 200GiB 10.8% 20
0GiB
ODEProb No Jac 1 42.0s 0.7% 42.0s 7.25GiB 0.4% 7.2
5GiB
ODE rhs Eval1 1 39.8s 0.7% 39.8s 4.62GiB 0.2% 4.6
2GiB
Create ODESys 1 19.1s 0.3% 19.1s 4.05GiB 0.2% 4.0
5GiB
Parse Network 1 15.6s 0.3% 15.6s 2.52GiB 0.1% 2.5
2GiB
DenseJac Eval1 1 2.44s 0.0% 2.44s 2.13GiB 0.1% 2.1
3GiB
SparseJac Eval2 1 103μs 0.0% 103μs 848B 0.0%
848B
ODE rhs Eval2 1 57.5μs 0.0% 57.5μs 752B 0.0%
752B
──────────────────────────────────────────────────────────────────────────
────
Picture of the solution
sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5) plot(sol, legend=false, fmt=:png)
